function r = chebev(a, b, c, m, x)
% Chebyshev evaluation: c is the array of T coefficients
% calculated by chebft. This function then evaluates the polynomial
% sum(k=1 to m) ck*tk(y) - c0/2, at the point y = (x-(b+a)/2)/((b-a)/2)
% this is the value of the function at x.
% From Press et al., p. 193. Coded Matlab P. Manis 21/Mar/99
%
r=NaN;
if((x-a)*(x-b) > 0)
   disp('x is not in range in chebev routine')
   return;
end
y = (2*x-a-b)/(b-a);
y2 = 2*y;
d = 0;
dd = 0;
for k=0:m-2
   j = m - k;
   sv = d;
   d = y2*d-dd+c(j);
   dd=sv;
end
r = y*d-dd+0.5*c(1);

